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ABSTRACT 


More  and  more  we  are  relying  on  satellites.  With  the  advancements  in  technology,  we 
will  soon  be  using  satellites  for  things  never  dreamed  of  before.  In  today’s  world,  some  of  these 
advancements  will  be  tied  to  saving  money  in  other  areas.  This  forms  the  foundation  of  many 
projects  looking  into  offboard  sensor  technology.  The  goal  is  to  save  money  by  using  space 
assets  to  replace  or  supplement  sensors  on  aircraft,  ships,  or  any  other  system. 

The  author  has  worked  on  several  projects  relating  to  offboard  sensors.  The  software  we 
are  presenting  here  has  been  a  valuable  tool.  This  software  package  is  written  for  Microsoft 
Excel.  It  enable  the  user  to  enter  orbital  information  and  viewing  angles  of  a  satellite.  The 
program  will  then  provide  several  plots.  These  plots  include  the  satellite  ground  track,  the 
coverage  region  of  the  satellite  over  a  time  interval,  and  the  instantaneous  field  of  view  for  any 
time.  Additionally  the  user  can  enter  a  square  target  region  of  any  size  with  its  center  at  a 
particular  latitude  and  longitude. 

With  this  software  package  the  user  will  be  able  to  look  at  various  orbits  and  viewing 
characteristics  and  easily  see  the  impact  of  the  choices  on  the  coverage  of  different  geographic 
regions.  This  is  very  helpful  for  the  initial  stages  of  designing  an  orbit  for  a  satellite  system.  It 
is  also  helpful  in  analyzing  the  coverage  of  existing  satellite’s  for  different  regions.  Finally,  this 
software  package  can  be  useful  for  students  just  entering  the  study  of  Astrodynamics. 
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1.  Background 


With  the  advancement  of  space  technology,  we  are  performing  many  functions  more 
efficiently  from  space.  The  potential  to  save  money  by  accomplishing  tasks  from  space  is 
leading  us  to  reconsider  the  way  we  conduct  business.  For  example,  we  may  soon  find  GPS 
replacing  the  thousands  of  radio  navigation  aids  (TACAN,  VOR,  NDB,  etc.)  as  the  primary 
means  of  navigation  for  aircraft.  With  the  current  emphasis  on  cutting  or  limiting  government 
spending,  new  projects  may  need  to  look  to  space  for  cost  effective  solutions. 

One  such  area  that  is  getting  a  lot  of  attention  these  days  is  the  idea  of  offboard  sensors. 
Think  of  the  number  of  aircraft  and  ships  worldwide  which  have  the  same  types  of  sensors.  We 
have  radar,  electronic  sensors,  and  photographic  capabilities.  Yet  we  are  accomplishing  all  of 
these  functions  from  space  as  well.  Thus  we  should  examine  the  feasibility  of  using  space  assets 
to  gather  the  information  and  relay  this  ,o  different  users. 

If  we  can  get  the  same  quality  of  information  from  space  to  the  user  in  a  timely  manner, 
there  is  the  potential  to  save  money.  Take  radar  as  an  example.  Most  aircraft  have  a  radar.  It  is 
heavy,  requires  a  lot  of  maintenance,  and  is  not  cheap.  If  we  are  able  to  use  space  assets  to 
provide  some  or  all  of  the  radar  information  to  the  aircraft,  we  may  be  able  to  use  a  smaller  and 
cheaper  radar  on  the  aircraft,  thus  saving  money. 

The  following  pages  contain  the  development  of  a  software  package  which  may  be 
helpful  in  the  initial  stages  of  examining  an  offboard  sensor  problem.  This  package  was 
developed  for  a  project  relating  to  the  concept  of  offboard  sensors. 


2.  Introduction 


While  working  on  a  project  relating  to  offboard  sensors,  the  need  for  software  to  help 
with  the  design  or  evaluation  became  necessary.  The  project  was  to  design  a  prototype  system 
which  could  provide  timely,  tactical  information  that  is  key  to  fighting  and  winning  isolated, 
short  duration  battles  anywhere  in  the  world  below  70  degrees  latitude.  Basically  we  were 
looking  for  a  satellite  which  could  gather  and  relay  information  to  a  mobile  user  in  a  tactical 
theater  environment.  The  idea  was  to  gather  or  relay  photo  reconnaissance,  electronic 
intelligence,  weather,  news,  and  force  location  to  friendly  forces  The  requirement  for  the 
satellite  to  support  a  specific  theater  (launch  on  demand),  perform  photo  reconnaissance,  and  act 
as  a  communications  relay  made  the  orbit  design  very  difficult. 

This  software  provides  the  means  to  analyze  different  satellite  orbits  and  viewing  angles 
versus  a  given  target  region.  Thus  if  we  have  a  battle  or  target  region,  we  can  visually  compare 
the  .satellite  ground  track  to  the  target  region.  Additionally,  we  can  use  different  elevation  angles 
(representing  how  much  of  the  earth  the  satellite  can  “see”  at  a  particular  time)  in  order  to 
examine  when  different  satellite  functions  can  be  performed  relative  to  the  target  region. 

For  example,  lets  assume  our  satellite  can  only  take  pictures  up  to  a  45  degree  elevation 
angle,  but  communications  are  possible  up  to  a  5  degree  elevation  angle.  The  software  will 
allow  us  to  see  which  passes  allow  communications,  which  allow  both  photos  and 
communications,  and  which  passes  will  have  no  coverage.  Additionally,  we  can  step  through 
minute  by  minute  to  get  a  snapshot  or  instantaneous  coverage  view  of  the  satellite  versus  the 


target  region.  This  software  is  very  useful  in  examining  orbits  and  elevation  angles  and  the 
impact  of  changes  on  the  coverage  of  a  given  target  region. 

2.  Astrodynamics 

Sections  3.1  through  3,3  are  a  review  of  the  basics  of  Astrodynamics.  For  the  reade*' 
faniiliar  with  the  two-body  equations  of  motion,  classical  orbital  elements,  the  Kepler  problem, 
and  the  effect  of  perturbations  on  the  orbital  elements,  you  can  may  skip  to  section  3.4. 

3,1  Equations  of  Motion 

Before  we  can  write  software  to  generate  satellite  orbits,  ground  tracks,  and  coverage,  we 
must  mathematically  describe  a  general  satellite  orbit.  The  following  equations  from 
Finidamentats  of  Astrodynamics  form  the  starting  point  for  describing  2-body  motion  (that  of  a 
satellite  about  a  large  body  such  as  the  Earth):' 


?  +  ((.tr'’)?  =  0 

(1) 

r  =p/(l  +  e  cos  v) 

(2) 

p  =  a(  1  -  e^) 

(3) 

Equation  (1)  is  the  2-body  differential  equation  which  describes  the  motion  of  a  satellite  around 
a  large  body.  This  equation  assumes  the  only  force  acting  on  the  satellite  is  gravity,  the  ma.ss  of 
the  central  body  (the  Earth  in  our  case)  is  much  greater  than  the  mass  of  the  satellite,  the  central 
body  is  a  uniform  sphere  (that  is  it  can  be  treated  as  a  point  mass),  and  there  are  only  two  bodies 
in  the  system  (the  Earth  and  the  satellite)  Equation  (2)  is  a  solution  to  equation  (1)  known  as  the 


polar  equation  of  a  conic. ^  The  parameter  |i  is  known  as  the  gravitational  parameto'-  and  is  the 
product  of  the  universal  gravitation  constant,  G,  and  the  mass  of  the  major  attracting  body,  the 
earth  in  our  case.^ 

All  orbits  (with  the  former  assumptions)  follow  the  path  of  a  conic  section  (ellipse, 
parabola,  or  hyperbola).  We  will  be  examine  only  elliptical  orbits  for  this  analysis.  Looking  at 
equations  (2)  and  (3)  the  parameters  “a”  and  “e”  represent  the  size  and  shape  of  the  elliptic  orbit. 
The  parameter  v  is  the  current  position  on  the  ellipse.  To  fully  define  a  satellite  orbit  we  can 
use  the  “Classical  Orbital  Elements”;  the  semi-major  axis  “a”,  the  eccentricity  “e”,  the 
inclination  “i”,  the  right  ascension  of  ascending  node  (or  longitude  of  ascending  node)  fl , 
argument  of  perigee  ai,  and  true  anomaly  v  to  define  an  orbit.'’  For  certain  orbits  Q,  co,orv 
may  be  undefined  and  are  replaced  with  different  parameters  (see  Space  Mission  Analysis  and 
Design  for  a  full  description  of  the  orbital  elements),’  Figure  1  below  shows  the  relationship 
between  a  satellite  orbit,  the  orbital  elements,  and  the  Geocentric-equatorial  coordinate  frame,  or 
UK  frame  (the  figure  was  modified  slightly  for  appearance):® 


3.2  Orbital  Elements  to  Position  Vector 


Our  initial  task  is  to  transform  the  orbital  elements  to  a  position  vector  in  the  UK  frame. 
Once  we  have  the  satellite’s  position  vector  in  the  UK  frame,  we  should  be  able  to  transform 
this  to  a  latitude  and  longitude.  The  first  step  in  this  process  is  to  define  the  orbit  in  the  Perifocal 
or  PQW  frame.  The  PQW  coordinate  system  as  it  relates  to  the  UK  frame  is  shown  in  Figure  2 
below:’ 

Figure  2  Relationship  betwee*’  PQW  and  UK 


The  frame  is  defined  with  P  in  the  direction  of  perigee  from  the  center  of  the  Earth,  Q  is  in  the 
orbit  plane  90  degrees  from  P  in  the  direction  of  satellite  motion,  and  W  is  perpendicular  to  the 
orbit  plane  to  form  the  third  axis  of  a  right  hand  system. 


ft 


The  following  equations  relate  the  orbital  elements  to  the  position  and  velocity  vectors  in 
the  PQW  frame.  ** 


r  =  r  cos  V  P  +  r  sin  V  Q 


V 


[-sinvP  +  (e+  cos  v)Q] 


(4) 

(5) 


Given  the  orbital  elements  ‘  a”,  “e”,  and  v  and  using  equations  (2)  through  (5)  we  can  find  the 
position  and  velocity  vectors  r  and  v  in  the  PQW  frame. 

Now  we  need  to  convert  the  position  and  velocity  vectors  from  the  PQW  to  the  UK 
frame.  Thus  we  must  perform  a  coordinate  system  transformation.  As  Figure  1  shows,  the 
orbital  plane  is  related  to  the  UK  frame  by  the  angles  “i”,  O,  and  w.  Using  these  three  angles, 
we  can  transform  any  vector  defined  in  the  PQW  frame  to  a  vector  defined  in  the  UK  frame 
(see  Appendix  A  for  the  coordinate  system  transformation).*^ 
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3.3  Orbit  Propagation 

3,3,1  Kepler  Problem  * 

So  far  we  can  define  a  satellite’s  orbit  through  the  use  of  the  orbital  elements  and 
transform  a  set  of  elements  to  a  position  vector  in  the  UK  frame.  Now  we  must  find  a  way  to  g 

propagate  the  satellite  through  its  orbit.  That  is,  we  must  be  able  to  find  the  satellite’s  position 
vector  at  any  time  after  the  initial  given  position. 

ft 


7 


ft 


9 


9 


9 


9 


9 


The  problem  of  finding  a  satellite’s  position  in  an  orbit  after  some  period  of  time  is 
known  as  the  Kepler  problem.  There  are  several  sources  which  go  into  great  detail  deriving  the 
solution  to  the  Kepler  problem.  Here  we  will  go  right  into  the  formulas  and  algorithms  for 
solving  this  problem.  The  following  is  a  summary  of  the  information  from  Fundamentals  of 
Astrodynamics  and  the  reader  can  refer  to  chapter  4  of  this  book  for  a  more  detailed  analysis.'" 

Out  of  the  six  orbital  elements,  only  v,  the  true  anomaly,  changes  with  time  (using  the 
two-body  spherical  earth  assumptions).  Thus  we  need  an  algorithm  which  will  give  us  the  true 
anomaly  of  the  satellite  after  some  time  period  Thus  if  we  are  given  a  value  for  v  at  T^O,  we 
want  to  find  the  new  v  at  other  time 

The  following  equations  are  necessary  for  solving  this  problem; 


„  e  +  cos  V 

cos  E  = - - 

1  +  e  cos  V 

M  =  E  -  e  sin  E 


M I  =  M  +  nAT 


cos  Vj 


e  -  cos  E| 
e  cos  E|  -  I 


(b) 

(7) 

(S) 

(b) 

(10) 


The  new  terms  are:  eccentric  anomaly,  E;  mean  anomaly,  M,  mean  motion,  n;  and  the  elapsed 
time,  AT.  We  will  assume  the  initial  orbital  elements  define  the  orbit  and  satellite  location  at  a 
particular  time  T=0.  In  equations  (6)  through  (10)  the  variables  M,  E,  v  are  for  T^O  and  the 
same  variables  with  a  subscript  1  are  for  a  particular  time  AT  relative  to  the  initial  condition 
time.  Equations  (6)  through  (10)  appear  to  be  straight  forward.  We  calculate  the  mean  motion 
and  the  eccentric  and  mean  anomalies  for  T=0  from  the  initial  orbital  elements.  Then  using 


8 


equation  (9)  we  get  the  mean  anomaly  for  the  time  of  interest,  Mj.  Once  we  calculate  E,,  we 
can  easily  find  v,.  After  comparing  Mj  to  V|  we  can  choose  the  correct  V]  (remember  the 
inverse  cosine  function  will  provide  two  answers,  the  true  anomaly  and  eccentric  anomaly  will 
both  be  between  0  and  180  or  between  180  end  360  degrees).  Thus,  we  can  find  V|  for  any  time. 

Of  course  the  calculation  of  the  eccentric  anomaly  from  mean  anomaly  is  not  as  P  'sy  as 
the  reverse.  Referring  to  equation  (6),  we  can  see  there  is  no  way  to  isolate  the  eccentric 
anomaly.  This  is  knov;n  as  a  transcendental  function.''  One  way  to  solve  for  eccentric  anomaly 
given  the  mean  anomaly  is  to  iterate.  That  is,  pick  a  value  for  eccentric  anomaly  and  see  what 
value  for  mean  anomaly  you  get.  We  iterate  until  the  difference  from  the  given  and  calculated 
mean  anomalies  is  sufficiently  small.  We  can  implement  the  following  algorithm  easily  on  a 
computer: 

1  Calculate  E,  M,  n,  then  Ml  based  on  AT 

2.  Set  E]  =  M|  as  an  initial  guess 

3.  En|,;w  =Mi  +  eSin(Ei) 

4.  Error  =  Absolute  value  of  (Eni.;w  -  E]) 

5.  Set  E]  =  En|.;w  ,  this  is  now  the  new  estimate  for  E] . 

6.  Perform  steps  3  through  5  until  the  error  is  less  than  some  appropriately  small 


7.  Use  equation  (10)  to  calculate  v,.  If  M,  is  in  quadrant  1  or  2  then  v,  is 
the  angle  from  the  inverse  cosine  function.  If  M,  is  in  quadrant  3  or  4  then  we  must 
take  360  minus  the  angle  we  get  from  the  inverse  cosine  function  to  get  the  proper  Vj, 

Thus,  if  we  define  a  set  of  orbital  elements  at  some  epoch  time  (T=0),  we  can  find  Vj  for  any 
time  using  the  algorithm  above.  However,  some  of  our  two-body  assumptions  may  cause  errors 
in  the  orbit  propagation.  To  improve  our  orbit  propagator,  we  will  try  to  account  for  some  of 
these  errors  by  looking  at  perturbations. 

3.3,2  Perturbations 

As  we  saw  in  section  3.3,  the  earth  is  not  a  perfect  sphere.  This  not  only  presents  a 
problem  for  finding  latitude,  it  impacts  the  motion  of  the  satellite.  For  a  spherical  Earth,  the 
force  due  to  gravity  would  only  depend  on  the  distance  of  the  satellite  from  the  center  of  the  UK 
frame  (which  is  at  the  center  of  the  “spherical  earth”).  Since  the  earth  is  oblate,  the  force  due  to 
gravity  will  vary  depending  on  where  the  satellite  is  even  if  the  distance  from  the  origin  of  the 
UK  frame  is  the  same.  Thus  the  orbit  will  have  a  torque  resulting  from  the  offset  center  of 
gravity. 

The  result  on  the  two-body  equations  of  motion  is  that  the  orbit  will  precess.  The  effect 
on  the  orbital  elements  is  that  and  co  will  change  over  time.  The  following  equations  describe 
the  change  in  these  elements  over  time  to  a  degree  accurate  enough  for  our  purposes:'^ 


=  -1.5  n  J2  (Rg  /a)^(cosiXl  e^)  ^ 

=  -2.06474xl0''^a“’^''^(cosi)(l-e^)“^  (11) 

TOj^  =  0,75  n  J2(Re  /  a)^(4  -  5sin^  i)(l  -  e^)“^ 

s  1.03237xl0''*a“^^^(4-5sin^il(l-e^)"^  (12) 

The  J2  term  accounts  for  the  earth  oblateness  and  R^  is  the  Earth’s  equatorial  radius  (the  same  as 
a^  in  section  3.4).  There  are  other  J„  terms,  known  as  zonal  coefficients,  which  account  for  the 
uneven  distribu'  u  11  of  mass  throughout  the  earth.  However,  the  Jj  term  is  responsible  for  a 
majority  of  the  perturbations  for  Earth  orbiting  satellites. 

Thus  with  the  information  we  have  so  far,  we  can  easily  write  code  which  will  propagate 
a  satellite  through  its  orbit.  We  start  with  a  set  of  orbital  elements  for  time  T=0.  For  any  time  of 
interest  we  can  calculate  the  new  orbital  elements  Q,  m,  and  v.  Next  we  transform  these  new 
elements  to  a  position  vector  in  the  UK  frame.  Now  we  Just  need  a  way  to  find  the  subsatellit  ^ 
latitude  and  longitude.  With  a  series  of  latitude  and  longitude  coordinates  we  can  form  the 
satellite  ground  track. 
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3.4  Position  Vector  to  Geodetic  Latitude  and  Longitude 

Now  that  we  can  find  the  position  vector  from  a  set  of  orbital  elements  at  any  time,  we 
must  find  a  way  to  convert  the  UK  position  vector  to  a  latitude  and  longitude  on  the  earth.  For  a 
spherical  earth,  the  calculation  of  a  latitude  and  longitude  from  the  position  vector  is  almost 
trivial.  Unfortunately  this  is  not  the  case.  The  earth  is  ('blate.  Here  we  will  assume  the  earth  is 
an  ellipse  with  a  polar  radius,  b^,  of  6356.755  km  and  an  equatorial  radius,  ae,  of  6378,14  km.'’ 
The  latitude  for  a  spherical  earth  is  known  as  geocentric  while  the  latitude  for  an  elliptic  earth  is 
known  as  geodetic.  Figure  3  below  shows  the  relationship  between  geodetic  and  geocentric 
latitude.''' 


Figure  3  Geodetic  versus  Geocentric  Latitude 


a=geodetic  latitude  p=geocentric  latitude 


The  “X”  component  of  the  satellite’s  position  in  Figure  3  is  formed  by  taking  the  resultant  of 
the  I  and  J  components  of  the  position  vector.  The  “Z”  component  is  the  same  as  the  K 
component  of  the  position  vector.  That  is: 


Thus  we  need  to  find  a  way  to  calculate  a  latitude  and  longitude  from  the  components  “X”  and 
“Z.” 

First  we  must  define  the  angles  a  and  p.  The  geocentric  latitude,  P,  is  the  “angle 
between  the  equatorial  plane  and  the  radius  from  the  geocenter,’’  and  geodetic  latitude,  a,  is  the 
“angle  between  the  equatorial  plane  and  the  normal  to  the  surface  of  the  ellipsoid.’’  We  can 
see  from  Figure  3  that  the  geocentric  latitude  is  simply  the  inverse  tangent  of  Z/X.  However, 
since  most  maps  use  geodetic  latitude,  it  would  be  useful  to  find  a  way  to  calculate  a  from  the  X 
and  Z  position  vector  components. 

Suppose  we  have  a  point  (X,  Z)  which  is  on  the  ellipse  and  the  current  position  of  the 
satellite  is  (Xo,  Zo).  If  the  slope  of  the  line  connecting  the  two  points  is  also  normal  to  the 
ellipse  at  (X,  Z),  thi  n  this  line  should  form  the  angle  a.  Pursuing  this  line  of  reasoning,  the 
following  equations  may  prove  useful: 


equation  of  ellipse 

X"  z"  , 
a.  b. 

(14) 

slope  of  line 

Z-Zo 

X-Xo 

(15) 

a  Z 

Slope  of  normal  to  point  on  ellipse  — ^ — 

b,X 

(16) 

Equation  (14)  defines  the  earth  ellipsoid.  Remembering  from  basic  calculus  that  the  derivative 
of  a  function  is  the  tangent  or  slope  function,  we  can  find  the  slope  of  the  tangent  to  the  ellipse  at 
any  point.  We  can  also  remember  the  negative  inverse  of  a  slope  is  the  slope  of  the  line 
perpendicular  or  normal  to  the  original  line.  Thus  if  dZ/dX  is  the  tangent  to  the  ellipse,  -dX/dZ 
is  the  normal  to  the  ellipse  at  any  (X,  Z)  which  is  on  the  ellipse.  By  differentiating  equation  (14) 
and  rearranging  terms  to  get  -dX/dZ  we  arrive  at  equation  (16). 

If  the  location  of  the  satellite  is  (Xo,  Zu),  the  solution  to  our  problem  will  be  a  point  (X, 
Z)  which  is  on  the  ellipse  and  allows  equations  (15)  and  (16)  to  be  equal.  This  sounds  like  an 
easy  problem  since  this  would  give  us  two  equations  with  two  unknowns.  However,  when  we 
solve  for  one  variable  and  then  substitute  into  equation  (14),  we  end  up  with  a  quartic 
polynomial.  See  Appendix  B  for  the  quartic  solution  to  this  problem. 


Thus,  we  need  a  different  approach.  Since  we  arc  looking  for  a  computer  solution,  an 
iteration  algorithm  seems  appropriate.  First,  we  isolate  the  problem  to  the  first  quadrant  in  the 
XZ  frame.  This  is  done  by  first  remembering  from  equation  ( 1 3)  that  Xo  must  be  positive. 


Since  the  sign  of  Zo  only  determines  whether  the  latitude  is  north  or  south,  wc  can  limit 


ourselves  to  the  first  quadrant.  After  we  calculate  the  value  of  a,  we  can  then  go  back  to  Zo  to 
determine  if  the  lati'  'de  is  north  or  south. 
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Now  to  devise  an  algorithm  to  solve  for  a.  Suppose  we  start  off  with  a  first  quadrant 
point  (X,  Z)  which  is  on  the  ellipse  and  a  satellite  position  (Xo,  Zo).  Plugging  these  values  into 
equations  (15)  and  (16)  will  provide  the  slope  of  the  line  from  the  point  to  the  satellite,  or 
“slope”,  and  the  slope  of  the  line  normal  to  the  ellipse,  or  “normal”,  at  (X,  Z).  If  the  slope  is 
greater  than  the  normal  then  our  choice  for  X  was  too  large.  If  the  slope  is  lower  than  the 
normal  our  choice  of  X  is  too  low.  We  can  iterate  by  increasing  or  decreasing  the  value  of  X 
until  the  difference  between  the  normal  and  the  slope  is  very  small.  The  portion  of  the  computer 
code  which  accomplishes  this  algorithm  is  in  Appendix  B. 

Now  with  the  difficult  part  complete,  we  move  to  the  much  easier  task  of  finding  the 
longitude.  Earlier  we  found  the  X  component  of  the  XZ  frame  position  vector  from  the  I  and  J 
components  of  the  UK  position  vector.  If  we  measure  the  angle  from  the  1  axis  eastward  along 
the  equator  to  the  X  axis,  we  get  an  angle  0  known  is  “local  sidereal  time,”  see  Figure  4  below.''’ 


^'igure  4  Local  Sidereal  Time  and  XZ  Plane 


Normally  we  measure  longitude  east  or  west  of  the  Greenwich  meridian.  Here  we  will  measure 
the  longitude  east  of  Greenwich  with  the  variable  A,],;.  The  angle  from  the  I  axis  to  Greenwich  is 


known  as  “Greenwich  sidereal  time,”  0g,  The  following  equations  and  Figure  5  below  show 
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how  all  these  variables  relate  to  each  other: 


0  — 

0  =  0  +tng,(t-to)  +  ?CE 


(H 

(18 


Figure  5  GST,  LST,  and  East  Longitude  Relationship 


The  parametei  100  is  the  angular  velocity  of  the  earth  and  0g„  is  the  Greenwich  sidereal  time,  or 
location  of  Greenwich,  at  a  particular  reference  time  t„.  While  we  are  not  using  the  SEZ  frame 
shown  in  Figure  5,  we  can  view  the  origin  of  the  SEZ  frame  as  the  subsatellite  point. 
Additionally  the  Z  axis  would  pass  directly  through  the  satellite’s  position. 

From  our  UK  position  vector  we  can  easily  calculate  0  by  taking  the  inverse  tangent  of 
rj/r,.  The  sign  of  the  J  component  will  give  use  the  proper  quadrant.  We  will  make  the 
assumption  for  software  that  Greenwich  is  at  the  1  axis  at  time  T=0.  In  other  words,  we  will 
initialize  the  program  with  0g„  at  0.  With  this  assumption,  we  can  solve  for  for  any  time  t 
using  equation  (18). 

This  assumption  really  makes  the  program  more  user  friendly.  We  initially  select  a 
longitude  of  ascending  node,  fl,  by  picking  a  longitude  which  is  appropriate  for  our  target 
region.  After  all,  our  interest  is  in  how  the  satellite  ground  track  relates  to  some  point  on  the 
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ground.  Once  we  pick  the  initial  orbital  elements  which  relate  in  some  way  to  our  target  region, 
the  rest  will  take  care  of  itself  We  account  for  the  rotation  of  the  earth  after  our  starting  point 
using  the  parameter  tEl0 .  Thus  we  can  now  calculate  the  latitude  and  longitude  of  the  sub¬ 
satellite  point  from  the  UK  position  vector. 

3.5  Satellite  Coverage 

Generating  a  ground  track  from  a  set  of  orbital  elements  may  sound  complex,  but  it  is  not 
that  difficult.  A  rough  sketch  of  a  ground  track  using  only  two-body  assumptions  is  a  fairly  easy 
task.  A  more  interesting  and  useful  problem  is  plotting  the  region  which  is  within  the  satellite’s 
field  of  view.  The  equations  developed  so  far  form  the  necessary  basis  for  allowing  us  to  solve 
this  problem.  Here  we  wish  to  solve  two  problems.  First,  we  need  to  plot  the  swath  ground 
track.  The  swath  ground  track  plot  shows  the  sea.  covered  by  the  satellite  as  it  moves  through 
its  orbit.  Our  second  problem  is  to  be  able  to  plot  the  instantaneous  field  of  view  (IFOV)  at  any 
time.  For  both  of  these  problems  out  goal  is  to  show  the  relationship  between  the  swath  ground 
track  and  IFOV  to  a  specific  ground  target  region, 

3.5.1  Swath  Ground  Track 

What  is  a  swath?  The  swath  is  “the  area  on  the  surface  of  the  Earth  around  the  ground 
track  that  the  satellite  can  observe  as  it  passes  overhead.’’”*  Since  a  picture  is  worth  a  thousand 
words,  we  will  plot  with  a  ground  track  and  a  swath  (the  region  between  the  outer  lines)  to  help 
illustrate.  In  Figure  6  below  is  a  plot  of  the  ground  track  and  swath  for  the  first  90  minutes  of  a 


satellite  orbit. 


Figure  6  Ground  Track  and  Swath  Coverage 
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We  can  see  from  this  plot  that  any  latitude  and  longitude  that  is  within  the  swath  will  be  visibl 
to  the  satellite  at  some  time  during  this  pass.  Now  we  must  develop  the  math  necessary  to 
generate  the  swath  plot. 


First,  we  must  define  some  angles  which  relate  the  satellite’s  view  of  the  earth  to  the 
swath.  Below  is  Figure  7  which  shows  the  geometry  and  angles  we  will  need  for  our 
calculations.''^ 


rhe  angle  X  is  the  earth  central  angle,  This  happens  to  be  half  the  swath  in  terms  of  earth  angle. 
The  angle  a  is  the  elevation  angle,  or  the  angle  above  the  horizon  the  satellite  is  to  an  observer 
on  the  earth.  This  will  be  one  of  the  input  parameters  for  the  software.  A  communications 
satellite  may  have  an  e  of  perhaps  5  to  10  degrees  while  a  photo  reconnaissance  satellite  might 
be  limited  to  30  to  45  degrees.  The  nadir  angle,  or  p  the  angle  off  the  subsatellite  point  which  is 
to  the  edge  of  the  field  of  view.  Finally  the  angle  p  is  the  angle  from  the  subsatellite  point  to  the 
horizon. 

Initially  we  will  look  at  the  angle  X.  From  this  angle  we  should  be  able  to  find  a  way  to 
generate  the  swath.  The  following  formulas  will  help  us  find  X  from  the  e  and  the  position 
vector  of  the  satellite.^" 


Re 

(19) 

sinp  - 

Re+  H 

sin  p  =  cos  e  sinp 

(20) 

X  -  90  -  p  -  e 

(21) 
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In  equation  (19),  Re  is  the  radius  of  the  Earth  and  H  is  the  satellite  altitude  above  the  Earth. 

Since  Re  +  H  is  the  same  as  the  magnitude  of  the  position  vector,  we  have  all  the  information 
necessary  to  calculate  k. 

Now  we  need  a  way  to  find  the  location  of  the  swath  for  a  given  satellite  position.  Figure 
8  below  shows  the  geometry  of  the  earth  central  angle,  X,  with  a  satellite  “flying  into  the  page.” 


Figure  8  Swath  Vector 


If  we  can  find  the  vector  y  in  Figure  8,  we  can  simply  add  y  to  r  to  get  the  swath  vector.  If  we 
subtract  y  from  r  we  would  get  the  swath  vector  on  the  other  side  as  well.  Since  the  velocity 

vector,  V  is  into  the  page,  we  can  form  the  vector  h  by  taking  the  cross  product  of  f  and  v . 
Then  we  can  use  the  equations  below  to  find  y  , 


y  =  |r|taiiA. 


h 


h  =  r*v 


(22) 

(24) 
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Using  equations  (22)  through  (24)  we  can  form  two  swath  vectors  by  adding  or  subtracting  y  to 
or  from  r.  Since  both  of  these  vectors  are  in  the  UK  frame,  we  can  find  the  latitude  and 
longitude  where  these  vectors  pass  through  the  earth  using  the  algorithms  from  section  3.4. 
Plotting  a  series  of  such  points  will  provide  us  a  swath  ground  track. 

3.5.2  Instantaneous  Field  of  View 

Not  only  do  we  want  a  swath  ground  track,  but  we  would  like  to  have  the  instantaneous 
field  of  view  for  any  particular  time.  The  swath  ground  track  will  allow  us  to  compare  a  target 
region  on  the  earth  with  satellite  coverage  over  time.  The  instantaneous  field  of  view  (IFOV) 
will  allow  us  to  analyze  the  satellite  coverage  for  any  time  and  compare  this  with  our  target 
region.  This  will  also  make  it  easy  to  visually  see  how  the  coverage  varies  for  a  particular  pass 
by  stepping  through  different  IFOVs. 

This  calculation  is  fairly  easy.  Since  we  are  dealing  with  a  specific  time,  we  will  have  a 
specific  satellite  position  vector.  From  the  methods  of  section  3.4,  we  can  find  the  subsatellite 
latitude  and  longitude  for  the  time  in  question.  The  following  equations  will  allow  us  to 
calculate  the  latitude  and  longitude  (5i,  L  |)  of  any  point  given  the  azimuth  (AZ),  earth  central 
angle  X,  and  the  subsatellite  latitude  and  longitude  (8,s,  Ls).^' 


AL  =  |L,  -L.,.| 

(25) 

cos(5.,  ')=  cosX.sin8j.  1- .sin  Xc.nsfi.(  cosAZ  ( 

5, V  180) 

(26) 

cosAL  =  (cosA  -  sin5„  sin6 , )  /  (cos5„  cos5  , ) 

(27) 

6.,.'=  90-  6.,. 

(28) 

Using  equations  (25)  through  (28)  and  varying  the  azimuth  angle  from  0  to  360  degrees,  we  can 
generated  a  series  of  points  which  are  at  the  edge  of  the  satellite’s  IFOV. 


4.  Software 


Now  that  we  have  all  the  equations  and  algorithms  from  section  3,  we  have  all  the  tools 
necessary  to  write  computer  code.  We  want  code  which  will  allow  us  to  generate  plots  of  the 
satellite  ground  track,  swath  ground  track,  IFOV,  and  even  the  orbit  relative  to  the  earth  as 
viewed  from  the  positive  W  direction  of  the  PQW  frame.  In  order  to  make  the  software 
available  to  many  users  and  possibly  students,  we  are  using  Visual  Basic  modules  for  Microsoft 
Excel.  Additionally,  the  subroutines  can  be  accessed  as  user-defined  functions.  Thus  the  user 
can  use  the  algorithms  to  solve  other  spreadsheet  type  problems.  Now  we  will  introduce  the 
software. 

41  Excel  Software 

We  are  using  Microsoft  Excel  mainly  for  the  ease  of  generating  user-defined  functions. 
While  the  user  may  find  this  software  slow  on  some  machines,  the  author  could  not  find  a  more 
widely  available  software  package  to  solve  this  problem.  Additionally,  the  ease  of  generating 
plots  with  Excel  and  importing  these  plots  to  documents  made  it  the  software  of  choice. 

Here  we  will  introduce  the  code  for  the  excel  file  PLOT.XLS.  All  the  code  for  PLOT  is 
in  Appendix  C.  Module  1  contains  all  the  number  crunching  code  and  algorithms.  Module  2 
contains  the  code  for  the  “bells  and  whistles,”  or  the  user  friendly  dialog  boxes.  All  of  the  plots 
are  made  by  calling  certain  functions  from  module  1  into  several  separate  spreadsheets.  We  then 
form  different  charts  by  using  the  data  from  the  spreadsheets.  The  modules  and  spreadsheets  for 
chart  generation  are  all  hidden  when  the  user  opens  the  PLOT.XLS  file. 
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PLOT  is  useful  for  short  time  spans  and  low  Earth  orbits  since  it  will  only  generate 
approximately  75  subsatellite  points.  For  higher  orbits,  or  if  the  user  desires  several  revolutions 
for  the  plots,  we  are  including  SLOWPLOT.XLS.  SLOWPLOT  is  the  same  program,  however 
it  will  plot  500  points.  As  the  name  implies,  it  is  slow.  It  will  take  rbout  one  minute  to  generate 
a  plot  from  a  set  of  orbital  elements  on  a  486-66  machine.  PLOT  will  only  require  about  10 
seconds  to  generate  plots  since  it  is  plotting  fewer  points. 

4.2  Software  Setup 

To  get  going,  simply  open  the  file  PLOT.XLS  or  SLOWPLOT  XLS  under  Excel  The 
user  will  find  five  chart  sheets  and  one  input  sheet  when  first  opening  the  file.  Once  the  usei 
opens  the  Excel  file,  it  is  very  important  to  select  the  manual  calculation  mode.  From  the 
toolbar  select  the  “Tools”  menu  bar  then  select  “Options.”  Click  on  the  “Calculation”  tab  and 
select  manual  calculation.  If  the  user  leaves  the  calculation  mode  in  automatic.  Excel  will  run 
though  all  the  calculations  every  time  the  user  changes  a  number.  Since  many  calculations  are 
necessary,  this  really  slows  down  the  program.  There  is  a  button  with  the  label  “Too  slow?” 
which  will  remind  the  user  what  steps  to  take  to  set  the  calculation  mode  to  manual.  This  button 
is  on  the  “Info”  sheet,  just  click  the  button  and  the  instructions  will  appear. 

Again,  once  the  file  is  open,  there  should  be  five  sheets  or  charts.  Each  has  a  title  on  the 
lower  tab.  The  sheets  and  their  contents  are: 

1 .  Gnd  Track  &  Inst  Cov  -  This  chart  shows  the  target  region,  the  satellite  ground 
track,  and  the  IFOV  for  any  specific  time. 


2.  Inst  Cov  Vs  Target  -  This  chart  shows  the  IFOV  for  any  specific  time  and  the 


target  region. 

3.  Swath  Plot  -  This  chart  shows  the  satellite  ground  track,  the  swath  track,  and 
the  target  region. 

4.  Gnd  track  -  This  chart  only  shows  the  ground  track  of  the  satellite. 

5.  Orbit  -  This  is  a  plot  showing  the  orbit  relative  to  the  earth.  It  shows  the 
satellite  position  at  T=0  and  the  satellite  position  for  any  specific  time. 

6.  Info  “  The  is  the  sheet  where  the  user  enters  the  data.  The  user  may  enter  the 
data  directly  into  the  spreadsheet  or  use  the  “buttons”  to  the  right  of  the  data  which  will  cause 
the  software  to  prompt  the  user  for  the  data. 

Sheets  !,  2,  and  5  have  a  “View  new  time  instant”  button  on  each  chart.  If  the  user  clicks 
this  button,  the  software  will  prompt  the  user  for  a  new  time  instant.  This  button  is  very  useful 
for  stepping  the  satellite  through  a  coverage  pass  of  a  region.  The  IFOV  for  different  time 
instants  can  provide  the  user  with  knowledge  of  how  the  coverage  changes  with  time  for  the 
pass. 

4.3  Input  Data 

Before  the  software  can  generate  any  plots,  the  user  must  enter  some  information.  The 
minimum  information  is  the  orbital  elements  for  T=0,  the  start  and  finish  times  relative  to  T=0 
for  the  plots,  and  an  instantaneous  plot  time  The  target  region  information  is  option;  1  If  the 
user  desires  a  target  region,  the  user  must  enter  the  latitude  and  longitude  for  the  center  of  a 
square  target  region.  Additionally,  the  user  must  enter  the  size  of  one  side  of  the  square  region 
in  kilometers. 
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With  respect  to  the  target  region,  there  is  one  other  optional  piece  of  information.  The 
user  can  select  the  number  of  orbits  between  target  region  plots.  This  value  is  set  to  0  if  the  user 
does  not  enter  a  number.  This  feature  will  plot  the  target  region  location  for  two  points  prior  to 
and  after  the  actual  target  region  location.  The  spacing  between  the  target  region  plots  is  the 
number  of  orbits  the  user  enters.  We  will  discuss  this  more  in  section  0  when  we  show  some 
sample  plots. 


In  Figure  9  below  is  a  sample  of  what  the  information  sheet  looks  like  with  actual  orbit 
and  target  region  information. 


Figure  9  Input  Information  Sample 
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0 
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95 
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60 
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1 

Calculated  Period 

100.876 

min 
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» 
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4.4  Sample  Plots 


Using  the  information  from  Figure  9,  we  will  now  show  what  the  some  of  the  plots  look 
like.  First,  we  will  start  with  the  chart,  the  “Gnd  I’rack  &  Inst  Cov”  in  Figure  10  below. 


Figure  10  Ground  Track  &  Instantaneous  Coverage  Plot 
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We  can  easily  see  from  this  plot  that  the  target  region  is  not  visible  to  the  satellite  at  7  60 
minutes,  However,  the  region  looks  as  if  it  was  visible  to  the  satellite  earlier. 

If  we  look  at  the  swath  plot  in  f  igure  1 1  below,  we  can  easily  see  the  target  region  is  in 


view  of  the  satellite  during  this  pass. 


Figure  1 1  Swath  Plot 
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In  I’act,  wc  can  sec  the  target  region  was  in  view  for  at  least  the  two  prior  passes  and  will  be  in 
view  the  next  two  passes.  I’he  target  regions  with  the  dashes  arc  where  the  target  region’s 
location  is  relative  to  the  satellite  orbit  on  other  passes.  Since  we  are  using  tine  orbit  spacing  for 
this  example,  we  can  see  the  location  ol'the  target  region  for  two  passes  prior  and  after  this  pass. 
The  spacing  takes  into  account  the  rotation  of  the  earth  and  nixie  regression  during  the  orbit. 

The  function  “Target”  in  Appendix  C  is  the  subroutine  which  perfonns  the  target  region 


calculations. 


One  nice  feature  of  the  program  is  the  ability  to  view  different  time  instants  easily.  We 
simply  click  the  “View  new  time  instant”  button  and  select  the  time  of  interest.  We  did  this  for 
T=20  minutes  in  the  “Inst  Cov  Vs  Target”  plot  in  Figure  12  below. 


Figure  12  IFOV  Vs  Target  Region  for  Specific  Time 
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We  can  easily  see  from  Figure  12  that  the  entire  target  region  is  in  view  of  the  satellite  at  T  20 
minutes. 


Another  plot  is  the  “Gnd  Track”  plot.  We  provide  ground  track  plot  in  case  the  user  only 
wants  the  ground  track  and  not  the  target  region,  swath,  or  II'OV  information.  Finally,  the  last 
plot  is  the  orbit  relative  to  the  earth.  Since  this  plot  shows  little  useful  information  for  circular 


orbits,  we  will  show  the  plot  of  a  Molniya  orbit.  In  Figure  13  below  is  the  plot  for  a  Molniya 


orbit  with  orbital  elements  of  a=26,600,  e=0.75,  i=63.4  deg,  oj=270  deg,  Q=100,  and  v=0.‘^ 


Figure  13  Orbit  Relative  to  Earth 
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This  is  one  orbit  that  PIOI'  does  not  provide  enough  points.  Sinee  the  satellite  is  Hying  very  fast 
near  perigee,  PLO  T  will  miss  much  of  the  ground  track.  Here  is  where  SLOWPLOT  can  be 
useful.  Since  SLOWPLOT  uses  more  points,  the  ground  track  will  not  have  the  largi  rap 
between  points  that  we  would  have  with  PLOT’.  Using  SLOWPLOT’  we  can  generate  the  Swath 
Plot  for  this  same  Molniya  orbit  (see  Figure  14  below). 
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Figure  14  MoIniyaPlotUsingSLOWPLOT.XLS 
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Now  that  we  now  the  tyfies  ol' plots  that  PI-OT  and  SLOWPLOT  can  generate,  we  should  move 
on  to  some  limitations. 


4.5  Limitations 

I'he  main  limitation  of  the  software  is  that  it  will  only  work  for  elliptical  orbits.  1  he 
program  will  not  plot  any  orbit  which  is  not  elliptical  (circular  orbits  are  obviously  okay  since 
they  are  just  a  special  case  of  the  ellipse).  Additionally,  the  program  is  only  as  good  as  its 
assumptions  In  summary.,  we  are  assuming  a  two  body  problem  with  an  oblate  Harth.  We  use 


I 


an  ellipsoid  model  for  the  Earth  with  a  polar  radius  of  6356.755  km  and  equatorial  radius  of 
6378. 14  km.  We  are  using  a  value  of  398600.5  km^/sec^  for  the  value  of  Also,  the  latitude 
and  longitude  calculations  rely  on  the  accuracy  of  the  spherical  earth  model.  Any  differences  for 
a  specific  location  from  the  ellipsoid  model  should  be  considered  as  necessary. 

Since  we  set  our  focus  on  satellite  orbits  which  are  “looking”  at  regions  of  the  Earth, 
orbits  above  geosynchronous  (GEO)  will  not  have  the  correct  perturbation  calculations.  Our 
assumptions  only  consider  the  perturbations  from  the  Earth’s  oblateness.  Thus  third  body 
effects,  such  as  the  earth  or  moon,  which  normally  become  noticeable  above  GEO,  are  not 
considered. 

These  limitations  really  have  no  impact  on  the  graphical  data  from  the  various  plots.  The 
errors  are  so  small,  and  the  plots  are  of  such  a  broad  scope,  the  impact  is  negligible,  We 
mention  these  for  the  case  where  the  user  decides  to  use  some  of  the  functions  for  other 
applications.  We  will  discuss  some  possible  applications  in  the  next  section, 


(§) 
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4.6  Applications 

There  are  many  possible  applications  where  this  software  can  be  useful.  Many  of  the 
functions  or  subroutines  in  the  software  have  been  used  to  solve  many  project  related  questions. 
This  software  is  a  great  as  an  aid  in  preliminary  design.  Different  concepts  for  an  orbit  can 
easily  be  examined  using  this  program.  Additionally,  this  software  can  be  very  useful  for 
students  learning  Astrodynamics.  Getting  a  visual  picture  of  how  the  orbit,  ground  track,  and 
satellite  coverage  changes  as  we  adjust  the  orbital  parameters  can  be  an  invaluable  tool. 
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Here  is  a  very  interesting  application.  Suppose  you  want  to  run  a  simulation  with  several 
satellites  in  orbit,  and  you  want  to  analyze  the  coverage  of  a  specific  region  over  time,  It  could 
be  possible  to  do  this  all  in  Excel  by  modifying  the  functions  and  adding  a  few  more.  But  lets 
say  you  want  to  use  a  package  similar  to  SIGMA  (Simulation  Graphical  Modeling  and  Analysis 
system)  to  generate  the  satellite  latitudes  and  longitudes  over  time.  We  did  this  with  six  lines  of 
code  in  one  node  by  using  a  variation  ofPLOT.XLS.  First  we  eliminate  the  perturbation  effects. 
Then  generate  a  data  file  with  one  period  of  latitudes  and  longitudes  and  the  time  for  each.  With 
this  data  file  you  can  generate  any  circular  orbit  with  the  same  inc  lination  and  semi-major  axis. 

It  is  a  simple  matter  to  calculate  node  regression  as  a  function  of  time.  With  a  time  correction  to 
account  for  the  true  anomaly,  and  a  longitude  correction  to  account  for  node  regression,  earth 
rotation,  and  different  Q's,  we  can  calculate  the  satellite  latitude  and  longitude  for  any  time. 

Thus  the  tools  are  here.  With  a  little  imagination  and  some  minor  changes,  they  can 
provide  solutions  to  many  problems.  Again,  the  main  goal  here  is  to  provide  some  help  with  the 
initial  designs  of  an  orbit.  Use  these  tools  to  get  you  in  the  ballpark,  then  use  the  more  powerful 
programs  for  the  refinement. 

5.  Conclusion 

Many  design  problems  require  finding  a  satellite  orbit  which  can  satisfy  different 
constraints.  The  software  packages  PLOT.XLS  and  SLOWPLOT.XLS  allow  us  to  visually 
examine  the  impact  of  various  orbit  parameters  on  the  satellite’s  coverage.  Such  a  tool  can  be 
invaluable  when  trying  to  balance  out  competing  parameters.  We  can  easily  see  how  changes  to 
the  design  will  impact  coverage  of  selected  regions  on  the  Earth.  Finally,  the  user  can  access  the 


various  functions  in  the  software  for  use  in  other  applications.  These  tools  have  proven  to  be 
very  helpful  in  analyzing  various  offboard  sensor  projects  in  the  past.  It  is  up  to  you  to  find  the 
uses  in  the  future. 
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Appendix  A  Coordinate  Transformation  PQW  to  UK 


Here  we  will  show  the  transformation  matrices  for  converting  a  vector  from  the  PQW 
frame  to  the  UK  frame.  As  Figure  2  from  section  3.2  shows,  the  PQW  frame  relates  to  the 
UK  frame  through  the  angle  Q,  i,  and  co.  Looking  at  Figure  2,  we  see  we  must  first  rotate  the 
PQW  frame  by  the  angle  -co  about  the  W  axis  (or  3  axis).  Then  we  rotate  the  PQW  frame  about 
the  P  axis  (or  1  axis)  by  the  angle  -i.  Finally  we  rotate  the  frame  one  more  time  about  the  W 
axis  by  an  angle  -Q.  This  will  result  in  the  two  frames  being  coincident.  These  rotation  are 
accomplished  using  rotation  matrices.  The  following  is  the  mathematical  process  for 
transforming  a  PQW  vector  to  an  UK  vector. 

‘10  O' 

0  cos(0)  sin(0)  =  rotation  about  1  axis  by  angle  0 

0  -sin(0)  cos(0) 

cos(P)  sin(P)  0" 

-sin(P)  cos(P)  0  ^  rotation  about  3  axis  by  angle  p 
0  0  1 

1  0  0  cos(-i)  sin(-i)  0  10  0  rp 

0  cos(-Q)  sin(-Q)  -sin(-i)  cos(-i)  0  0  ccis(-^)  sin(-^)  rg 
0  -  sin(-r2)  cos(-r2)  0  0  1  ^  '  sin(-Tn)  cos(-'m)_ 


Appendix  B  Geodetic  Latitude  Solutions 


Quartic  Solution 

The  following  is  the  derivation  of  the  quartic  solution  for  the  problem  of  converting  a 
satellite  position  to  geodetic  latitude.  First  we  convert  the  UK  position  vector  to  the  XZ  plane: 

Xo  =  Zo  =  tk 

Then  we  find  a  point  (X,  Z)  on  the  Earth  ellipsoid  such  that  the  slope  of  the  line  form  (Xo,  Zo) 
to  (X,  Z)  is  the  same  as  the  slope  of  the  line  normal  to  the  ellipse  at  (X,  Z). 


Equation  for  the  Earth  ellipsoid: 


,4.2  i).2  ■  ’ 


Equation  for  the  slope  of  line  normal  to  the  ellipse  (-dX/dZ) : 


Ae^Z 

Be'^X 


Slope  of  line  from  point  on  ellipse  to  satellite: 


Z-Zo 

X-Xo 


Set  the  two  slope  equations  equal  to  each  other  and  isolate  X: 


Z-Zo  _Ae^Z 
Be^X 


XoAe^Z 


{Be^- Ae'^)z-ZoBe'^ 


Now  square  X  piuy  inio  ellipse  equation: 


{.{Be^'-Ae'^)z-ZoBe-y 


Xo^Ae^Z^  I  ^1 

{jyBt?- -  Ae^^  Z  -  ZoBe^^ 


Solve  for  Z: 


{Be^-2  Be^Ae'^+Ae^)  +  i2Ae‘^  ZoBe"^  -2Be^  ZoBe^) 
+  (zo;  '  ■^Xo^Ae'^Be'^-Be'^Ae^-Be^  +  lBe^Ae'^) 

^{iBe^  ZoBe^  -  2  Be^Ae^  ZoBe'^)  Z-Be^  ZoBe"^  =  0 

The  correct  Z  will  be  a  real  solution  to  the  quartic  which  has  the  same  sign  as  Zo.  After  we  find 

Z  it  is  a  simple  matter  of  substitution  to  find  X.  Finally,  the  geodetic  latitude  is  the  inverse 

tangent  of  (Z-Zo)/(X-Xo), 


The  following  is  the  Visual  Basic  code  used  in  the  modules  for  the  Excel  program  to 
solve  this  problem  using  iteration 


'TO  FIND  GEODETIC  LATITUDE  FROM  r  IN  ijk  FRAME 
Function  LA1(R1,  RJ,  RK) 

Dim  r,  X,  Z,  SLP,  NOR,  INC  As  Single 

STP  -  0 

CHN  =  0 

AE- 6378. 14 

BE  =  6356.785 

r  =  mag(RI,  RJ,  RK) 

INC  =  3# 

'CHECK  FOR  EAT  CLOSE  TO  EQUA'I  OR  OR  POLE 
'There  is  little  ditTerence  between  geodetic  and  geocentric  latitude 
'close  to  these  extremes.  Here  we  will  assume  they  are  equal  when 
'within  5  deg  of  the  equator  or  pole. 

If  (Abs(ASlNR(RK  /  r)  180  /  PlO)  <  5)  Or  (Abs(ASINR(RK  /  r)  180  /  Pl())  >  85)  I  hen 
LAT-ASlNR(RK/r)'^  180/P1() 

Else 

'We  will  start  with  an  initial  guess  for  X  then  iterate  to  .solve 
'start  with  X  for  geocentric  ca'culation  as  the  initial  guess 
X  =  AB*Sqr{RI''2  i  RJ''2)/r 
Z  -  Sqr(BE  ^  2  (1  -  X  ^  2  /  AE  ^  2)) 

SLP  (Z  -  Abs{RK))  /  (X  -  Sqr(RI  2  +  RJ  2)) 

NOR  -(Z'^AE"'2/(X*BE''2)) 

If  SLP  >  NOR  Then  CHN  =  I  'X  is  too  large  so  set  tlag  to  decrease  X 
Do 

If  AbsCSLP  -  NOR)  >  0.0001  Then 
'slope  too  small  so  increase  X 
If  ((SLP  <  NOR)  And  CHN  -  0)  Then  X  X  +  INC 


iS 


■M 


'slo  ie  too  small  but  we  decreased  X  last  iteration 
If  ((SLP  <  NOR)  And  CHN  =  1 )  Then 
CHN  =  0 
INC  =  INC  /  3# 

X  =  X  +  INC 
End  If 

'slope  too  large  and  we  increased  X  last  iteration 
If  {(SLP  >  NOR)  And  CHN  -=  0)  Then 
CHN=  1 
INC  =  INC  /  3# 

X  =  X  -  INC 
End  If 

'slope  too  large  so  decrease  X 
If  ((SLP  >  NOR)  And  CHN  =  1 )  Then  X  =  X  -  INC 
Else 

'error  is  small  enough  so  calculate  latitude 
LAT  Atn(NOR)*  180/P1() 

'check  for  north{H )  or  south(-)  latitude 
If  (RK  <  0)  Then  LAT  =  -LAT 
STP  =  1 
End  If 

'X  should  never  be  larger  than  radius  of  Earth 
If  X  >  AE  Then  X  -  AE 
Z  =  Sqr(BE  2  '*  ( 1  -  X  ''  ?  /  AE  2)) 

SLP  =  (Z  -  Abs(RK))  /  (X  -  Sqr{Rl  ^  2  +  RJ  ^  2)) 

NOR  =  Z'*  AE''2/(X  '*  BE^2) 

Loop  While  STP  -  0 
End  If 

End  Function 

Contparision  of  Analytic  with  Iterative  Solution 

When  wc  call  this  function  a  position  vector  of  r=(7000  1,  7000  J,  7000  K.),  we  get  a 
latitude  of  35.361 13,  This  is  using  an  error  between  the  normal  and  the  slope  of  0.0001 .  On  the 
next  page  we  have  the  analytical  solution  to  the  quartic  equation  using  Maple  V.  Using  Maple  V 
we  get  a  latitude  of  35,35972729  degrees.  1  his  represents  a  0,004%  error.  Since  the  errors  from 
our  other  a.ssumptions  will  probably  introduce  more  error  than  this,  we  should  be  safe  at  this 
accuracy. 


(S) 


i 


i 


I 


i 


I 


» 


» 


B-J 


Maple  V  Solution  For  Specific  Example 


Position  vector  UK  frame  (7000,7000,7000) 
Ae:=6378. 1 4:Be:=6366.755:Xo:=9899.495;Zo:=7000: 

Ae  :=  6378.14 
Be  :=  6356.755 
Xo  :=  9899.495 


Zo  :=  7000 

ellipse:=X''2  /Ae''2  +  Z'^2  /  Be^2=1  ;slope_contraint:=  (Z-Zo)/(X-Xo)-(Ae''2/ 
Be''2)*ZyX: 


ellipse 


Ae~  Be^ 


.  .  Z-  Zo  Ae^  Z 

slope  contramt  :=  — — -  -  — ^ — 

X-Xo  Be-X 

s:=solve<{ellipse,slope_contraint},{X,Z}); 

,V  :={z=-.  1564955393  10^,  .4396945462  10^}, 

{Z  =  3670.495837,  5207.446490 }, 

{ Z  =  -3653.2435 1 2,  X  =  -52 1 9.64570 1  } 

lat;=arctan((Zo-3670.495837)/(Xo-5207.446490))*180/3.141592654: 

lat  :=  35.35972729 


Appendix  C  Visual  Basic  Code 
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Module  1:  Calculation  Code  For  PLOT.XLS 


'to  find  the  magnitude  of  a  vector 
Function  mag(a,  B,  C) 
mag  =  Sqr(a  '*  a  +  B  B  +  C  C) 

Bnd  Function 

'to  find  the  cross  product  of  two  vectors 
Function  CROSS(a,  B,  C,  D,  e,  F) 

CROSS  Array(B  F  -  C  e,  C  D  -  a  F,  a  c  -  B  D) 

End  Function 

» 

'to  find  the  dot  product  of  two  vectors 
Function  DOT(a,  B,  C,  D,  o,  F) 

DOT  =  a  *  D  ^  B  ♦  e  +  C  '♦  F 
End  Function 

'to  find  the  angle  between  two  vectors 
l-'unction  ANGLE(a,  B,  C,  D,  e,  I') 

ANGLE  180  /  IM()  *  ACOSR(DOT(a,  B,  C,  D,  e,  F)  /  (mag(a,  B,  C)  *  mag{D,  e,  F))) 

End  Function 

'to  find  the  inverse  cos  in  radians  FOR  THIS  MODULE  ONLY 
Function  ACOSR(X) 
lf(-X'*X+  l)<-0Thcn 
ACOSR  0 

Else  j 

ACOSR  -  Atn(-X  /  Sqr(-X  X  +  I ))  +  Pl()  /  2 
End  If 

End  Function 

'to  find  the  inverse  SIN  in  radians  FOR  1  HIS  MODULE  ONLY  ^ 

Function  ASINR(X) 

If  (-X  X  +  I )  <-  0  Then 
ASINR=  PI{)/2 
Else 

ASINR  Atn(X  /  Sqr(-X  X  r  I )) 

End  If 

End  Function 
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TO  FIND  RPQW  FROM  ORBITAL  ELEMENTS 
Function  RPQW(a,  e,  truean) 

Dim  ANG,  r  As  Single 

ANG  -  truean  *P1()/  180 

r“a*(l  -e''2)/(l  +  e*Cos(ANG)) 

RPQW  -  Array(r  *  Cos(ANG),  r  *  Sin(ANG),  0) 

End  Function 

TO  FIND  VPQW  FROM  ORBLI  AL  ELEMENTS 
Function  VPQW(a,  e,  truean) 

Dim  ANG,  P,  Cl  As  Single 
ANG  =  truean  *  Pl()  /  180 
P-  a*(l  -c'^2) 

Cl  -(398600.5  /P)'' 0.5 

VPQW  -  Array(-Cl  *  Sin(ANG),  Cl  *  (e  +  Cos(ANG)),  0) 
End  Function 


'PERFORM  COORDINATE  TRANSFORMATION  ANGLES  IW  DEGREES 
'AXIS  1  ROTATION 
Function  RO  l(a,B,C,ANG) 

Dim  X  As  Single 
X-ANC^PIO/ISO// 

ROTl  =  Array(a,  B  ♦  Cos(X)  t-  C  *  Sin{X),  -B  *  Sin(X)  t  C  *  Cos(X)) 

End  Function 

'AXIS  2  ROTATION 
Function  ROT2(a,  B,  C,  ANG) 

Dim  Y  As  Single 
Y  -  ANG  '^PI()/  180/( 

ROT2  -  Array(a  ♦  Cos(Y)  -  C  *  Sin(Y),  B,  a  ♦  Sin(Y;  i  C  Cos(Y)) 

End  Function 

'AXIS  3  ROTA  TION 
Function  ROT3(a,  B,  C,  ANG) 

Dim  Z  As  Single 
Z  =  ANG  PIQ/  180# 

ROT3  =  Array(a  Cos(Z)  +  R  "'  Sin(Z),  -a  *  Sin(Z)  +  B  "■  Co.s(Z),  C) 

End  Function 

Function  PI() 

PI -3. 141 59265358979 
End  Function 
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'PERFORM  COORDINA  TE  TRANSFORMATION  ANGLES  IN  DEGREES 
'  THE  FIRST  ROTATION  IS  ABOUT  THE  3  AXIS  BY  A  DEGREES 
'  THE  SECOND  ROTATION  ABOUT  THE  I  AXIS  BY  B  DEGREES 
'  THE  THIRD  ROl'ATION  ABOUT  THE  3  AXIS  BY  C  DEGREES 

Function  ROT313(RP,  RQ.  RW,  a,  B,  C) 
temp  -  ROT3(RP,  RQ,  RW,  a) 

TEMP2  ■  ROTI{temp(0),  temp(l),  temp(2),  B) 

ROT3 1 3  ROT3(TEMP2(0),  TEMP2(  I ),  TEMP2(2).  C) 

End  Function 


"FO  FIND  THE  VECTOR  DEFINING  THE  SWATH  PERPENDICULAR  I'O  GROUND 
TRACK 

Function  SWATH(RI,  RJ,  RK,  VI,  VJ,  VK,  EL) 

Dim  elrad,  r,  nadir,  lambda,  MAGY,  MAGH,  C  As  Single 
cirad  EL'»  PI()/ 180 
r  “  mag(FiL  RJ,  RK) 

nadir ASlNR(Cos(elrad)  6378. 145  /  r) 
lambda  ^  Pl()  /  2  -  elrad  -  nadir 
MAGY  -  r  3'an( lambda) 

H  =  CROSS(Rl,  RJ,  RK,  VI,  VJ,  VK) 

MAGH  --  mag(H(0),  H(l),  H(2)) 

C- MAGY /MAGH 

SWATH  --  ArraylC  H(0),  C  H(  1 ),  C  *  H(2)) 

I/nd  Function 

'This  function  will  take  any  angle  in  degrees  and  reduce  it  to  an  angle 
'between  0  and  360  degrees 
Function  MODULUS(numer,  denom) 

Dim  REALVALUE  As  Single 
'perform  integer  division  first 
intvalue  numer  \  denom 
REALVALUE  -  numer  /  denom 

MODULI  >S  (Abs( REALVALUE)  -  Abs( intvalue))  '*  denom 
If  REALVALUE  <  0  Then  MODULUS  -  -MODULUS 
End  Function 

'This  function  will  solve  the  transcendental  equation  from  the  Kepler 
problem,  That  is,  given  mean  anomaly  and  eccentricity, 

'find  the  eccentric  anomaly 
Function  ECCNEW(M1,  e) 

Dim  ERR,  El  As  Single 


lil  =M1 
Do 

ECCNEW  =  M1  I  e*  Sin(El) 

ERR  =  Abs(ECCNEW  -El) 

El  =ECCNEW 
Loop  Until  (ERR  <0.0001) 

If  Ml  >  PI()  And  ECCNEW  <  PI()  Then  ECCNEW  =  2  *  Pl{)  -  ECCNEW 
End  Function 


'This  will  solve  the  Kepler  problem 

'THIS  WILL  RETURN  THE  VALUE  OF  THE  NEW  TRUE  ANAMOLY  AFTER  TIME  'IN 
SECONDS 

Function  KEPLER(a,  e,  TNOM,  time) 

Dim  V,  ENOM,  N,  M,  MNEW,  VNEW,  ENEW  As  Single 
v=  TNOM '*P1()/  180# 

ENOM  -  ACOSR((c  i  Cos(v))  /  (I  +  e  '»  Cos(v))) 

If  V  >  PI()  And  ENOM  <  PI()  Then  ENOM  =  2  Pl()  -  ENOM 
N  -Sqr(398600,5/a^3) 

M  -  ENOM  -  e  ♦  Sin(ENOM) 

MNEW  -  M  +  N  '*  time  'TIME  MUST  BE  IN  SECONDS 
MNEW  =  MODULUS(MNEW,  2  PI()) 

ENEW  -  CCCNEW(MNEW,  e) 

VNEW  1 80  /  Pl()  ♦  ACOSR((c  -  Cos(ENEW))  /  (e  •  Cos(ENEW)  -  1 )) 

If  (ENEW  >  Pl()  Or  (ENEW  >  -Pl()  And  ENEW  <  0))  And  VNEW  <  180#  Then 
KEPLER  =  360#  -  VNEW 
Else 

KEPLER  - VNEW 
End  If 

End  Function 


'this  will  convert  orbital  elements  to  position  and  velocity  vector  in 

'the  UK  frame  for  given  time  T  in  minutes 

Function  rvattime(semi,  e,  i,  OMEGA,  ARGPiilR,  truean,  T) 

Dim  TRUENEW,  CONST  1,  WNEW,  ASCENNEW  As  Single 
TRUENEW  KEPLER(semi,  e,  truean,  T  ♦  60) 

PQWR  =  RPQW(semi,  e,  TRUENEW) 

POWV  -  VPQW(semi,  e,  TRUENEW) 

CONSTl  -  (1  -(e)  ■''2) ''(-2) 


'THIS 
!fe  = 


WILL  CORKECI 
0  Then 


’  FOR  PERTUBAl'iONS  T  IS  IN  MINUTES 


WNEW  -  0 


Else 

WNEW  -  ARGPER  +  71692361 1 1 1#  T  (semi)  ^  (-3,5)  (4  -  5  *  (Sin(i  "■  PI()/  180))^ 
2)  ♦  CONS  T  1 
End  If 


} 


I 


C-4 


ASCENNEW  =  OMEGA  -  14338472220#  *  T  ♦  (semi)  (-3.5)  ♦  Cos(i  *  PI()  / 180)  * 
CONST  1 

TJow  change  to  UK  frame.do  first  rotation  with  -WNEW,  argument  of  perigee 
'do  the  second  rottion  with  -inclination 

'do  the  third  rotation  with  ASCENNEW,  the  longitude  of  acsending  node 
rijk  =  ROT313(PQWR(0),  PQWR(l),  POWR(2),  -WNEW,  -i,  -ASCENNEW) 
vijk  -  ROT313(PQWV(0),  POWV(l),  PQWV(2),  -WNEW.  -i,  -ASCENNEW) 
rvattime  =  Array(rijk(0),  rijk(  I ),  rijk(2),  vijk(O),  vijk(  1 ),  vijk(2)) 

End  Function 

'this  will  generate  the  Latitude  and  Longitude  for  the  subsatellite  point 
'and  the  swath  for  the  time  interval  T1  to  TF  (minutes)  for  the  given  elevation 
'angle  elev.  Time  is  considered  0  at  the  given  orbital  elements. 

Function  ORBlT(semi,  e,  i,  OMEGA,  ARGPER,  truean,  Tl,  TF,  elev) 

Dim  TEST(74,  6),  TRUENEW,  CONSTl,  WNEW,  ASCENNEW,  INC,  T  As  Single 
lNC-(TF-TI)/73# 

T '  -  Tl 

For  N  =  0  To  73 

temp  =  rvattime(semi,  c,  i,  OMtiGA,  ARGPER,  truean,  T) 
rijk  Array(temp(0),  temp(  1 ),  temp(2)) 
vijk  =  Array(temp(3).  temp(4),  temp(5)) 

'Now  find  the  two  vectors  defining  the  swath  perpendicular  to  the  ground  track 
Y  -  SWATI  l(rijk(0),  rijk(  I ),  rijk(2),  vijk(O),  vijk(  I ),  vijk(2),  elev) 

SWI  -  AiTay(rijk(0)  i  Y(0),  rijk(l)  f  Y(l),rijk(2)  t  Y(2)) 

SW2  =  Array(rijk(0)  -  Y(0),  rijk(  1 )  -  Y(  I ),  rijk(2)  -  Y(2)) 

'Now  calculate  the  subsatcllite  point  (SSP)  latitude  and  longitude 
'  for  this  time  T 

'GREENWICH  MliRlDl AN  IS  INITIALIZED  AT  THE  1  AXIS  AT  TIME  i) 

TEST(N,  0)  -  T 

TEST(N,  1 )  "=  LON(rijk(0),  rijk(  I ),  rijk(2),  T) 

TES  f(N,  2)  =  LAT(rijk(0),  rijk(  I ),  rijk(2)) 

'FIND  THE  LONG  AND  LAT  OF  FIRST  SWATH  VECTOR 
TEST(N,  3)  -  LON(SW1(0),  SWI ( 1 ),  SWI (2),  T) 

TEST(N,  4)  -  LAT(SW1(0),  SWl(l),  SWI (2)) 

'FIND  THE  LONG  AND  L.A'f  OF  SliCOND  SWA  TH  VECTOR 
TEST(N,  5)  -  LON(SW2(0),  SW2(  I).  SW2(2),  T) 

TEST(N,  6)  =  LAT(SW2(0),  SW2(  I ),  SW2(2)) 

T  =  T  +  INC 
Next  N 

ORBIT  -  TEST 
End  Function 


'THIS  WILL  FIND  THE  INSTAN  TANEOUS  FIEI.D  OF  VIEW  FOR  THE  SA  TELLITE 
'A  T  A  PARTICULAR  TIME  ( MINI  I  TES) 


Function  instant(semi,  e,  i,  OMEGA,  ARGPER,  truean,  elev,  T) 

Dim  fovloc(73,  2),  slat,  slon,  elrad,  r,  nadir,  lambda  As  Single 
temp  =  rvattime(semi,  e,  i,  OMEGA,  ARGPER,  truean,  T) 
ri  jk  =  AiTay(temp{0),  temp(l),  temp(2)) 
slon  =  LON(rijk(0),  rijk(l),  rijk(2),  T) 
slat  =  EAT(rijk(0),  rijk(l),  rijk(2))  *  Pl()/ 180# 
elrad  =  elev  *  Pl()  /  180 
r  =  mag(rijk{0),  rijk{l ),  rijk(2)) 
nadir  =  ASINR(Cos(elrad)  *  6378. 145  /  r) 
lambda  ="  Pl( )  /  2  -  elrad  -  nadir 

'first  find  the  start  with  the  LCjN  AND  LAT  of  the  SSP  for  this  time  instant 
fovtoc(0,  0)  =  slon 
fovtoc(0,  l)  =  slat*  180/PI() 

■CAirULATE  THE  LON  and  LAT  AROUND  THE  SSP  AT  1  HE  DESIRED  SWATH 
X=  1 

For  AZ  -  0  To  360  Step  5 

tlat  =“  PI()  12  -  ACOSR(Cos(lambda)  *  Sin(slat)  ^  Sin(lambda)  *  Cos(slat)  *  Cos(AZ  *  Pl( )  / 
180#)) 

DELTALON  -  ACOSR((Cos(lambda)  -  Sin(slat)  Sin(tlat))  /  (Cos(s!at)  Cos(tlat)))  *180/ 

PIO 

fovloc(X,  l)-tlat*  180/P1() 

If  AZ<=  180  Then 
fovloc(X,  0)  =  slon  +  DEL, TALON 
Else 

fovloc(X,  0)  -  slon  -  DELTALON 
End  If 
X  =  X^  1 
Next  AZ 
instant  =-  fovloc 
End  Function 

'THIS  WILL  generate:  I'HE  1  ARGET  GRIDS  FOR  THE;  DESIRED  TIMES 
'  I'he  input  parameter  "num"  is  the  number  of  satellite  orbits 
'between  target  grid  plots.  Semi  is  the  semi-major  axis  of  the  orbit 
'which  is  necessary  to  determine  then  period  of  then  satellite 
Function  target! tlat,  tion,  size,  semi,  num,  e,  i) 

Dim  COOR(24,  I ),  RAD,  PER,  DEGCH,  LUP,  LLOW,  CENTER  As  Single 
RAD  (size  /  6378. 145)*  90 /Pl() 

PER  -  0,0001658669  *  semi  (1.5) 

'calculate  node  regression  in  degrees  per  orbit 

nodereg  = -23849472.41  *  semi  ^  (-2)  *  Cos(i  *  PI() /  180)*  (I  -  e  2) (-2) 

'now  calculate  how  far  the  earth  rotates  per  satellite  orbit 
DEGCH  -  PER  *  num  *  0.2506844773 
LUP  =  tlat  +  RAD 
LLOW  tlat  -  RAD 


x  =  o 

For  T  =  -2  To  2 

'calculate  the  latitude  of  the  square  grid 
COOR(X,  1)  =  LUP 
COOR(X+  1,  l)  =  LLOW 
COOR{X  +  2,  l)=LLOW 
COOR(X  +  3,  1)  =  LUP 
COOR(X  +  4,  1)  =  LUP 

'calculate  the  center  target  longitude  is  for  this  satellite  pass 
CENTER  -  tlon  +  T  (DEGCH  +  nodereg) 

'get  the  center  longitude  between  0  and  360  degrees 
CENTER  =  MODULUS(CENTER,  360) 

If  CENTER  <  0  Then  CENTER  =  CENTER  +  360 
COOR(X,  0)  =  CENTER  -  RAD 
COOR(X  +1,0)-  CENTER  -  RAD 
COOR(X  +  2,  0)  -  CENTER  i  RAD 
COOR(X  +  3,0)  =  CENTER  +  RAD 
COOR(X  I  4,  0)  =  CENTER  -  RAD 
X  =  X  +  5 
Next  T 

target  =  COOR 
End  Function 

'this  function  calculates  the  earth  centrral  angle  between  two 

'sets  of  latitude  and  longitude 

Function  MAPANGLE(slat,  slon,  tlat,  tlon) 

SET  =  slat '*P1()/  180 
Sin  =  slon  '*  Pl()/ 180 
TLT  -  tlaf^PlQ/  180 
TLN  =  tlon  Pl()/  180 

MAP  ANGLE  =  ACOSR(Sin(SLT)  *  Sin(TLT)  +  Cos(SLT)  *  Cos(TLT) '»  Cos(Abs(Sln  - 
TLN)))*  180/PI() 

End  Function 


I 


ft 


ft 


ft 


ft 


"fO  FIND  GEODETIC  LATITUDE  FROM  r  IN  ijk  FRAME 
Function  LAT(R1,  RJ,  RK) 

Dim  r,  X,  Z,  SLP,  NOR,  INC  As  Single 

STP  =  0 

CHN  =  0 

AE-6378.14 

BE  =  6356.785 

r  -  mag(RL  RJ,  RK) 

INC  -  2ti 


ft 


•  •  • 


•  •  •  • 


•  • 
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'CHECK  FOR  LAT  CLOSE  TO  EQUATOR  OR  POLE 
'There  is  little  difference  between  geodetic  and  geocentric  latitude 
'close  to  these  extremes.  Here  we  will  assume  they  are  equal  when 
'within  5  deg  of  the  equator  or  pole. 

If  (Abs(ASlNR(RK  /  r)  ♦  180  /  PI())  <  5)  Or  (Abs(ASlNR(RK  /  r)  *  1 80  /  PI{))  >  85)  Then 
LAT  =  ASINR(RK  /  r)  180  /  PI() 

Else 

'We  will  start,  with  an  initial  guess  for  X  then  iterate  to  solve 
'start  with  X  for  geocentric  calculation  as  the  initial  guess 
X  =  AE  *  Sqr(RI  2  +  RJ  ^  2)  /  r 
Z-Sqr(BE'^2'*‘(l  -  X '' 2  /  AE  ^  2)) 

SLP  -  (Z  -  Abs(RK))  /  (X  -  Sqr(Rl  ^  2  +  RJ  a  2)) 

NOR  =  (Z  *  AE  ^  2  /  (X  *  BE  2)) 

If  SLP  >  NOR  Then  CHN  =  1  'X  is  too  large  so  set  tlag  to  decrease  X 
Do 

If  Abs(SLP  -  NOR)  >  0.0001  Then 
'slope  too  small  so  increase  X 
If  ((SLP  <  NOR)  And  CHN  -  0)  Then  X  -  X  +  INC 
'slope  too  small  but  we  decreased  X  last  iteration 
If  ((Sl,P  <  NOR)  And  CHN  ■=  I )  Then 
CHN  -  O 
INC  -  INC  /  3# 

X-X  )  INC 
End  If 

'slope  loo  large  and  we  increased  X  last  iteration 
If  ((SLP  >  NOR)  And  CHN  =-  0)  I'hcn 
CHN  -  1 
INC  -■  INC  /  2ii 
X  -  X  -  INC 
End  If 

'slope  too  large  so  decrease  X 
If  ((SLP  >  NOR)  And  CHN  I )  Then  X  X  -  INC 
Else 

'error  is  small  enough  so  calculate  latitude 
LAT  -  Atn(NOR)*  180/PI() 

'cheek  for  north(  ^ )  or  south(-)  latitude 
If(RK  <0)Then  LAT  -LAT 
STP  1 

I*' 

l:,iiu  1 1 

'X  should  never  be  larger  than  radius  of  Earth 
IfX  >  AEThenX  =  AE 
Z  =  Sqr(BH  ^  2  *  ( I  -  X  ^  2  /  AE  ''  2)) 

SLP  =  (Z  -  Abs(RK))  /  (X  -  Sqr(RI  ^  2  +  RJ  ''  2)) 

NOR  -  AE'^2/{X*  BE  2) 

Loop  While  S  I  P -0 


End  If 

End  Function 

TO  FIND  GEODETIC  LONGITUDE  FROM  r  IN  UK  FRAME 
Function  LON(RI,  RJ,  RK,  TIM) 

IfRl  =  OThen 
If  RJ  >  0  Then  LON  =  90 
If  RJ  <  0  Then  LON  =  270 
Else 

LON  =  Atn(RJ  /  R!)  *  1 80  /  PI() 

If  RI  <  0  Then  LON  =  LON  +  1 80 
If  R1  >  0  And  LON  <  0  Then  LON  =  360  +  LON 
End  If 

CORLON  =  LON  -  0.2506844773  ♦  TIM  '  CORRECT  FOR  ROTATION  OF  THE  EARTH 
If  CORLON  <  0  Then 
Do 

CORLON  =  CORLON  +  360 
Loop  Until  CORLON  >0 
End  If 

LON  =  CORLON 
End  Function 

Module  2:  Dialog  Box  Code  for  PLOT.XLS 


Private  a,  c,  i,  om,  w,  v,  TI,  TF,  tins,  LT,  LNG,  size,  numorb,  elev 
Sub  INFOO 
default  -  0 

msgl  ■  "Welcome  to  the  orbit  plotter.  " 

msg2  "Greenwich  is  initialized  at  0  DEG  at  Time  T^O.  " 

msg3  “  "Now  we  will  get  some  information  about  the  satellite  orbit.  " 

msg4  =-  "IS  THE  ORBIT  CIRCULAR?" 

msg  =  msgl  ^  msg2  h  msg3  +  msg4 

Style  vbYesNo 

I:  response  ■■  MsgBox(rnsg,  Style,  "ORBIT  INFORMATION") 

If  response  vbYes  Then 

ALT  InputBoxC'ENTER  THE  ALTITUDE  (KM)", ,  default) 
a- ALl' t  6378.145 
e  0 

perigee  ^  a  *  ( 1  -  e) 

If  perigee  <  6378. 145  Then 

msgl  -  "This  orbit  has  a  perigee  which  will  impact  the  earth.  " 
msg2  ^  "  Retry  or  quit,  this  orbit  cannot  be  plotted" 
response  ^  MsgBox(msgl  &  msg2,  5,  "WARNING!") 

If  response  vbRetry  Then 
Gofo  1 


» 


» 


» 


i  • 


» 


C-9 


•  • 


m| 


i  k 


»  1 
w  I 


1 


response  =  MsgBox(msgl  &  msg2,  5,  "WARNING!") 

If  response  =  vbRetry  Then 
GoTo  1 
Else 
GoTo  2 
End  If 
End  If 

i  =  InputBoxC'ENTER  THE  INCLINATION". ,  default) 

If  i  ^  0  Then 
om  =  0 
w  =  0 

V  =  InputBoxC'ENTER  THE  TRUE  LONGITUDE  AT  T=0", ,  default) 

Else 

om  =  InputBoxC'ENTER  THE  LONGITUDE  OF  ASCENDING  NODE", ,  default) 
w  =  0 

V  =  InputBoxC'ENTER  THE  ARGUMENT  OF  LATITUDE  AT  T=0",  .  default) 
End  If 

Else 

a  =  InputBoxC'ENTER  THE  SEMI-MAJOR  AXIS  (KM)", ,  default) 
e  =  InputBoxC'ENTER  THE  ECCENTRICITY  ", ,  default) 
perigee  =  a  *  ( 1  -  e) 

If  perigee  <=  6378.145  Then 

msgl  =  "This  orbit  has  a  perigee  which  will  impact  the  earth.  " 
msg2  =  "  Retry  or  quit,  this  orbit  cannot  be  plotted" 
response  =  MsgBox(msgl  &  msg2,  5,  "WARNING!") 

If  response  =  vbRetry  Then 
GoTo  1 
Else 
GoTo  2 
End  If 
End  If 

i  =  InputBoxC'ENTER  THE  INCL.INATION", .  default) 

Ifi  =  0Then 
om  =  0 

w  =  InputB.'xC'ENTER  THE  LONGITUDE  OF  PERIAPSIS",  ,  default) 

V  =  InputBoxC'ENTER  THE  TRUE  ANOMALY  AT  T=0". ,  default) 

Else 

om  =  InputBoxC'ENTER  THE  LONGITUDE  OF  ASCENDING  NODE",  ,  default) 
w  =  InputBoxC'ENTER  THE  ARGUMENT  OF  PERIAPSIS",  ,  default) 

V  =  InputBoxC'ENTER  THE  TRUE  ANAMOLY  AT  T=0", ,  default) 

End  If 

End  If 

period  =  lnt(0.000 1658669  *  a '^(1.5)) 

msg9  =  "Your  selected  orbit  has  a  peril  f  " 

msglO  -  "Would  you  like  to  plot  only  oiit,  orbit  beginning  at  T=0?" 


•  •  • 


C-10 


msg  =  msg9  &  period  &  "  minutes.  "  &  msglO 
response  =  MsgBox(msg,  Style,  "PLOT  INFORMATION") 

If  response  =  vbYes  Then 
TI  =  0 
TF  =  period 
Else 

msgl  =  "The  times  asked  for  here  are  the  number  of  minutes  " 
msg2  =  "before  T=0  (negative  values)  or  after  T=0  (positive  values),  " 

TI  =  InputBox(msgl  &  msg2  &  "Enter  the  start  time  for  the  plot.",  ,  default) 
TF  =  InputBox("Enter  the  stop  time  for  the  plot  (in  minutes).", ,  default) 

End  If 

msgl  =  "One  plot  shows  a  snaphot  of  the  satellite  field  of  view  for  a  given  time. 

tins  =  InputBox(msgl  &  "Enter  the  time  (min)  for  the  snapshot  plot.",  ,  default) 

msgl  =  "Enter  the  minimum  elevation  angle  (deg)  desired  " 

msg2  =  "  (the  minimum  angle  from  the  earth  horizon  to  the  satellite).  " 

elev  -  lnputBox(msgl  &  msg2,  ,  default) 

LT  =  0 
LNG  =  0 
size  =  0 
numorb  =  1 

Worksheets("Info  Sheet"), Range("c5  ")  =  a 

Worksheets("Info  Sheet"), Range("c6  ")  =  e 

Worksheets("lnfo  Sheet"), Range("c7  ")  =  i 

WorksheetsC'Info  Sheet"). Range("c8  ")  =  om 

Worksheets("Info  Sheet"), Range("c9  ")  =  w 

WorksheetsC'Info  Sheet"). Range("clO  ")  =  v 

WorksheetsC'Info  Sheet"). Range("cl  1  ")  =  elev 

WorksheetsC'Info  Sheet").Range("cl3  ")  =  TI 

WorksheetsC'Info  Sheet" ),Rangc("c  14  ")  =  TF 

WorksheetsC'Info  Sheet").Range("cl5  ")  -  tins 

WorksheetsC'Info  Sheet").Range("cl8  ")  -  LT 

WorksheetsC'Info  Sheet"). Range("c  19  ")  =  LNG 

WorksheetsC'Info  Sheet"). Range("c20  ")  =  size 

WorksheetsC'Info  Sheet"). Range("c21  ")  =  numorb 

msgl  =  "Orbit  info  loaded.  You  can  now  enter  target  info  if  desired  " 

msg2  =  "or  click  on  PLOT  NOW  to  see  plots.  Click  OK  to  continue." 

r  =  MsgBox(msgl  ti  msg2,  0) 

2: 

End  Sub 


Sub  tgtINFOO 
Style  ^  vbYesNo 

msgl  =  "We  can  pul  a  square  grid  on  the  plots  representing  " 
msg2  =  "a  region  of  interest.  Simply  enter  the  latitude  and  longitude  " 


msg3  =  "for  tlie  center  of  the  region  and  the  desired  length  of " 

msg4  =  "the  grid  in  kilometers.  Click  YES  to  continue  or  NO  to  skip  this  feature." 

response  =  MsgBox(msgl  &  msg2  &  msg3  &  msg4.  Style,  "TARGET  INFORMATION") 

If  response  =  vbYes  Then 

LT  =  InputBox("Enter  the  region  Latitude  (north  is  +).") 

LNG  ^  InputBoxf  "Enter  the  region  Longitude  measured  East  of  Greenwich.") 

size  =  InputBoxC'Enter  the  length  of  the  square  region  (km).") 

nisgl  =  "We  can  also  plot  the  location  of  the  target  region  relative  to  " 

msg2  =  "the  ground  track  for  previous  and  subsequent  orbits.  " 

msg3  ^  "  Click  Yes  to  use  and  No  to  skip  this  feature." 

response  =  MsgBox(msgl  &  msg2  &  msg3.  Style,  "TARGET  INFORMATION") 

If  response  =  vbYes  Then 

nuniorb  =  InputBoxC'Enter  the  nui  ’er  of  orbits  between  relative  target  region  plots") 
Else 

numorb  =  0 
End  If 

Worksheets("Info  Sheet").Range("cl  8  ")  =  LT 
WorksheetsC'lnfo  Sheet"). Range("cld  ")  =  LNG 
Worksheets("Info  Sheet"). Range("c20  ")  -  size 
WorksheetsC'lnfo  Sheet"'  Range("c21  ")  =  numorb 
msgl  =  "Target  info  loaded.  " 

msg2  -  "  When  ready  to  generate  plots  click  on  Plot  Now.  " 
r  =  MsgBox(msgl  &  msg2,  0) 

End  If 
End  Sub 

Sub  INSTO 

tins  “WorksheetsC'lnfo  Sheet").Range("cl  5  ") 

WorksheetsC'lnfo  Sheet"). Range("c  15  ")  =  lnputBox("",  "Time  for  the  snapshot  plot  (min) ", 
tins) 

WorksheetsC'lnfo  Sheet"). Calculate 
WorksheetsC'SHEETl  (2)"). Calculate 
WorksheetsC'SHEET  1  (5)"). Calculate 
End  Sub 

Sub  PLOTNOWO 

a  -  WorksheetsC'lnfo  Sheet"), Range("c5  ") 
e  =  WorksheetsC'lnfo  Sheet"). Range("c6  ") 

Debug. Print  e 

If  a  <"  0  Or  e  <  0  Or  e  >=  1  Then 
msgl  =  "This  orbit  is  not  elliptical.  " 
msg2  ^  "  Sorry,  this  program  can  only  plot  elliptic  orbits." 
response  ^  MsgBox(msgl  &  msg2,  0,  "WARNING!") 

GoTo  3 


perigee  =  a  (1  -  e) 

If  perigee  <=  6378. 145  Then 

msgl  =  "This  orbit  h  perigee  which  will  impact  the  earth.  " 

msg2  "  Check  the  orbital  elements  and  make  the  necessary  corrections." 

response  =  MsgBox(msgl  &  msg2,  0,  "WARNING!") 

GoTo  3 
End  If 

Application. Calculate 
SheetsC'Swath  Plot"). Select 
3: 

End  Sub 
Sub  welcomeO 

msgl  =  "This  program  will  take  orbital  element  information" 

rnsg2  =  "  and  make  several  plots.  The  satellite  ground  track,  swath,  " 

msg3  =  "instantaneous  coverage  and  position  in  orbit  around  the  Earth  " 

msg4  =  "for  any  entered  time  can  be  displayed." 

r  =  MsgBox(msgl  &  msg2  &  msg3  &  msg4,  0) 

msgl  "=  "Additionally  we  can  plot  a  target  region.  We  can  then  easily" 

msg2  -  "  compare  the  ground  track  and  satellite  position  at  a  given  time" 

msg3  =  "  to  the  target  region." 

r  =  MsgBox(msgl  &  msg2  &  msg3,  0) 

msgl  =  "For  new  users  the  best  way  to  begin  is  to  click  on  the  ORBIT" 
msg2  =  "  INTERVIEW  button  and  answer  all  the  questions.  Then  click  on  " 
msg3  =  "  the  ADD  TARGET  REGION  button  if  you  want  to  enter  a  target  region. 
msg4  •  "Finally  click  on  PLOT  NOW  to  generate  the  plots." 
r ""  MsgBox(msg  1  &  msg2  &  msg3  &  msg4.  0) 

msgl  =  "  fhe  information  summary  table  contains  all  the  parameters  you  enter." 
msg2  =  "  You  can  change  any  value  directly  in  the  table  then  click  the" 
msg3  =  "  PLOT  NOW  button  to  generate  the  new  plots.  " 
r  =  MsgBox(msgl  &  msg2  &  msg3,  0) 

msgl  =  "Additionally  some  charts  have  a  VIEW  NEW  TIME  INSTANT  button." 
msg2  =  "  Clicking  this  button  moves  the  satellite  and  satellite  view  " 
msg3  =  "  to  the  new  time.  7'his  is  very  useful  in  comparing  the  satellite  " 
msgl  =  "  location  to  the  target  region  at  ditTerent  time  instants  " 
r  =  MsgBox(msgl  &  msg2  &  msg3  &  insg4,  0) 

msgl  =  "Depending  on  the  speed  of  your  machine,  it  may  take  a  minute" 
msg2  =  "  or  so  to  generate  the  plots.  Be  patient,  a  chart  will  appear " 
msg3  =  "after  all  calculations  are  complete.  There  v/il!  be  five  different " 
msg4  =  "charts,  just  click  on  the  chart  you  wish  to  see." 
r  =  M,sgBox(msgl  &  msg2  &  msg3  &  msg4,  0) 

End  Sub 

Sub  slow() 

msgl  '  "If  you  are  experiencing  long  delays  after  entering  data" 


msg2  =  "  (except  for  the  PLOT  NOW  button)  there  could  be  " 

msg3  =  "  an  easy  fix.  From  the  TOOLS  menu  select  OPTIONS.  Select 

msg4  =  "  Calculation.  Make  sure  the  calculation  mode  is  set  to  Manual. 

r  =  MsgBox(msgl  &  msg2  &  msg3  &  msg4  &  msgS,  0) 

msgl  =  "If  this  doesn't  help,  you  just  have  a  slow  processor." 

msg2  =  "  There  are  numerous  calculations  being  performed,  be  patient. 

r  =  MsgBox(msgl  &  msg2,  0) 

End  Sub 


